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In an effort to build a stronger microscopic foundation for radiation damage models in gallium ar- 
senide (GaAs), the electronic properties of radiation-induced damage clusters are studied with atom- 
istic simulations. Molecular dynamics simulations are used to access the time and length scales re- 
quired for direct simulation of a collision cascade, and density functional theory simulations are used 
to calculate the electronic properties of isolated damaged clusters that are extracted from these cas- 
cades. To study the physical properties of clusters, we analyze the statistics of a randomly-generated 
ensemble of damage clusters because no single cluster adequately represents this class of defects. The 
electronic properties of damage clusters are accurately described by a classical model of the electrical 
charging of a semiconducting sphere embedded in an uniform dielectric. The effective band gap of 
the cluster depends on the degree of internal structural damage, and the gap closes to form a metal 
in the high-damage limit. We estimate the Fermi level of this metallic state, which corresponds to 
high-energy amorphous GaAs, to be 0.46 ±0.07 eV above the valence band edge of crystalline GaAs. 
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[. INTRODUCTION 

Modeling the effect of radiation damage on an electronic 
circuit is a multiscale effort. Circuit simulations 1 rely on com- 
pact damage models^, which are constructed from device sim- 
ulations of single circuit elements 3 . In turn, electronic device 
simulation relies on defect models 4 that include the drift, dif- 
fusion, creation, and annihilation of defects with appreciable 
concentrations and non-negligible electronic activity. Identi- 
fication and characterization of such defects requires atom- 
istic simulations. These simulations are further divided into 
multiple scales and categorized by the choice of microscopic 
degrees of freedom. Binary collision approximation (BCA) 
simulations 5 only track atoms that deviate from the ideal crys- 
tal structure and are only limited by the total amount of crystal 
damage simulated rather than simulation volume. Classical 
molecular dynamics (MD) simulations^ track all atoms in a 
simulation volume, which restricts the simulation size while 
more accurately accounting for interatomic interactions. Den- 
sity functional theory (DFT) simulations 7 include a quantum 
description of electrons and enable the simulation of charge 
fluctuations, but with an asymptotic cost that scales cubicly 
with simulation volume. Fidelity of atomistic simulations is 
particularly important because microscopic errors can propa- 
gate through the simulation scales into macroscopic errors. 

We consider the case of gallium arsenide (GaAs), where 
modeling is necessary because of a relative dearth of experi- 
mental characterization of defects and radiation damage com- 
pared to silicon. The central assumption of radiation dam- 
age modeling at the device scale is that damage can be re- 
duced to a time-dependent spatial distribution of point de- 
fects and their binary chemistry with electrons, holes, and 
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other defects. In accordance, recent DFT studies have char- 
acterized the electronic properties of all likely point defects 
in GaAs^l The point-defect ansatz is consistent with the 
BCA of high energy collision cascades that consist primar- 
ily of well-separated Frenkel pairs, but it breaks down at the 
end of a collision cascade when a particle loses sufficient en- 
ergy that the mean distance between collisions approaches the 
interatomic separation (the "end of range"). MD simulations 
of this regime are consistent with the BCA in accounting for 
populations of point defects, but they also contain damage 
clusters with a comparable number of atoms that cannot be 
decomposed into a collection of isolated point defects^. 

The end of range in an MD simulation is the spatial re- 
gion surrounding the terminus of a single collision cascade 
event. Standard experiments do not resolve a single cascade 
event and the end of range instead refers to the average depth 
at which cascades terminate. The end-of-range layer in ion 
implantation experiments contains the highest density of ma- 
terial damage and nanoscale clustering of defects is observed 
regularlv 11 . The damage to GaAs caused by ion implantation 
and neutron irradiation produces broad anomalous features 
(the 'U' and 'L' bands) in deep-level transient spectroscopy 
(DLTS) 12 , which is an experimental probe used to character- 
ize the energy levels of point defects. Such features are absent 
in pristine and electron-irradiated GaAs. Experimental char- 
acterization of end-of-range damage remains an open problem 
and will benefit from microscopic insights derived from the- 
ory and simulation. 

The structural and electronic properties of GaAs damage 
clusters have not been studied to assess their importance in 
radiation damage models relative to point defects. For point 
defects in GaAs, there are a manageable number of plausible 
structures that can be enumerated and simulated exhaustively. 
In contrast, damage clusters can contain a large number of 
atoms and span a very large configuration space. No single 
structure would be adequately representative or predictive of 
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the properties of damage clusters. Instead, a study of their 
properties requires a statistical analysis of many cluster con- 
figurations. 

In this paper, we use a combination of MD and DFT sim- 
ulations to model the electronic properties of damage clus- 
ters in GaAs and strengthen the microscopic foundation upon 
which future radiation damage models will be built. First, we 
collect a sample set of damage clusters generated from MD 
simulations of collision cascades. Because of the discrepancy 
in computationally tractable simulation volumes between MD 
and DFT, the sampling is biased towards smaller clusters that 
are accessible to DFT calculations. We supplement this set 
with a larger sample set of artificially generated idealized 
damage clusters that are uniformly distributed over a range of 
cluster sizes. This artificial set is large enough to generate reli- 
able statistics, which enables us to carefully assess finite-size 
effects and extrapolate models into a regime of cluster size 
that is not directly accessible with DFT. To complement the 
extrapolation, large damage clusters are studied indirectly by 
generating bulk amorphous GaAs structures that are assumed 
to be representative of the cluster interior. The atomistic sim- 
ulation methodology is reviewed in Section[II] Several models 
of the electronic properties of defect clusters are proposed in 
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ported in Section IV and their implications for device-level 
radiation modeling are discussed in Section[V] 



II. COMPUTATIONAL METHODOLOGY 

The MD simulations are performed using the LammpjP 
code. They employ a Pettifor bond-order potential for GaAs 
developed by Murdick and co-workers^ that has been mod- 
ified at small interatomic separations to match the Ziegler- 
Biersack-Littmark (ZBL) universal repulsive potential^. The 
basic methodology for simulating displacement cascades is 
well-established and has been reviewed by Averback and de 
la Rubia 15 . In order to simulate a neutron-induced recoil, an 
atom is selected at random and given an impulse in a random 
direction with the desired energy. The system is evolved for a 
period of 10 ps using a variable time step which ensures that 
no atom moves more than 1 pm in a given step. The size of the 
simulation cell depends on the recoil energy to be simulated. 
The simulations of 1 keV recoils contain 216,000 atoms, while 
the simulations of 50 keV recoils contain 13,824,000 a toms. 
More details of these simulations are described elsewhe re 16 ! 17 ! 

Regions of GaAs that cannot be decomposed into point 
defects embedded in crystalline GaAs are categorized as 
amorphous regions. They are identified in the MD simulations 
using a generalization of an approach used earlier by Foiles 10 
to analyze displacement cascades in silicon. The criteria used 
to identify the amorphous regions is the presence of 5- and 
7-atom rings. A ring in this context is a path along nearest- 
neighbor atoms that returns to its original site. In the case of 
an ideal diamond structure, the smallest rings contain 6 or 8 
atoms along the path. To define an amorphous region, all rings 
of 7 atoms or less are identified. The system is divided into 
cells that correspond to the cubic cell of the diamond struc- 



ture, and the number of 5- and 7-atom rings centered in each 
cubic cell is determined. Cells that contain more than one 
5-atom ring or more than two 7-atom rings are classified as 
amorphous regions. 

The DFT simulation volume used in this study is a 4 x 4 x 4 
supercell of the conventional cubic unit cell of GaAs (a 11.2 
nm 3 cube) containing 256 Ga and 256 As atoms. Damage 
clusters are required to be contained within this volume with 
approximately 1 nm of crystalline GaAs as a buffer between 
the damage cluster and its periodic image. Amorphous re- 
gions that satisfy this criterion are visually identified and ex- 
tracted from the MD simulations. Ten damage clusters are 
generated using this procedure and then relaxed with MD to 
their metastable ground states. In three cases, the extracted 
subsystem is further annealed with MD for 100 ns at 300 K. 

Two bulk amorphous GaAs structures are also creating us- 
ing MD simulation, starting from a 4 x 4 x 4 cubic supercell of 
the ideal crystal. The crystal is heated to 2500 K and isother- 
mally evolved for 1 ns, then the temperature is reduced to 300 
K over 15 ns and isothermally evolved for 10 ns. This process 
produces a large number of Ga-Ga and As-As nearest neigh- 
bors, which are energetically unfavorable relative to Ga-As 
nearest neighbors. A second, lower-energy amorphous struc- 
ture is created from the first structure by applying a Monte 
Carlo procedure to pairwise exchange atoms in order to max- 
imize the number of Ga-As nearest neighbors. This structure 
is heated to 800 K, cooled to 300 K over 5 ns and isothermally 
evolved for 10 ns. 

Artificial damage clusters are generated using a simple 
procedure to simulate high-energy damage. We again begin 
the process with 4x4x4 cubic supercells of the ideal GaAs 
crystal. Within a sphere of a given radius that is arbitrarily 
centered on an As atom, all atoms on lattice sites and all va- 
cancies on interstitial sites are randomly assigned to new sites, 
lattice or interstitial, within the sphere. 5 different radii are 
chosen - 8, 10, 12, 14, and 16 bohr (1 bohr « 0.0529 nm) - 
and 25 clusters are generated for each of these radii, for a total 
of 125 artificial clusters. 

DFT calculations of the damage clusters are performed 
with the SeqQuest code 18 . The details of these calculations 
are guided by previous studies of point defects in GaAgSI The 
local density approximation (LDA) is used for electron ex- 
change and correlation 19 . Calculations using PBE 20 for sim- 
ple intrinsic defects do not exhibit any physically meaning- 
ful differences^! Norm-conserving pseudopotentials without 
semicore rf-orbitals but with nonlinear core corrections are 
used for both Ga and As atoms^H and a double-zeta-plus- 
polarization quality atomic orbital basis set is used to repre- 
sent crystal orbitals. The Brillouin zone is sampled only at the 
r point for the 4x4x4 supercells, sufficient for the accuracy 
needed in the current study. Gaussian smearing of 0.03 eV 
is applied to orbital occupations to aid in convergence of the 
self-consistent field cycle in the presence of multiple energy 
levels near the Fermi level. The lattice constant is fixed at the 
LDA value for GaAs, 0.560 nm. All damage clusters are re- 
laxed to the DFT local energy minimum using an accelerated 
steepest descent update of atomic positions to a tolerance of 
0. 1 eV/nm. Vertical charge transition energies are calculated 
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by constraining the structure of the charged damage clusters 
to the neutral geometries. Adiabatic charge transition energies 
are calculated by relaxing the structure of each charge state to 
its local energy minimum starting from the relaxed geometry 
of the previous charge state with one fewer unit of charge. 

DFT calculations of charged damage clusters use the lo- 
cal moment countercharge (LMCC) method 22 . The charge 
density in the unit cell is partitioned into a Gaussian approx- 
imation of the defect charge plus a charge-neutral remainder. 
The electrostatic potential within the unit cell is constructed 
from the periodic charge neutral remainder and the nonperi- 
odic Gaussian defect charge. The total energy is corrected 
with an embedding energy to align the charged defect to the 
perfect crystal potential and a bulk screening energy to ac- 
count for polarization outside the unit cell. A supercell proce- 
dure is used to account for bulk screening energy outside the 
volume of the unit cell, which was calibrated to converge to 
the bulk dilute limit in GaAs with a 1.6 bohr skin depth 8 . 

DFT calculations of the two bulk amorphous GaAs struc- 
tures are performed with the Vasp code 7 . The details of these 
calculations are slightly different from the SeqQuest sim- 
ulations. The LDA is still used to relax the structure, but 
the lattice constant is kept at the experimental value of 0.565 
nm. Default projector augmented wave (PAW) pseudopoten- 
tials are used for both Ga and As, without e?-orbitals in the 
valence (consistent with the SeqQuest calculations). A dif- 
ferent model of electron exchange and correlation, the Heyd- 
Scuseria-Ernzerhof (HSE06) modeffl, is used for electronic 
density of states (DOS) calculations. The resulting HSE06 
band gap of crystalline GaAs is 1.26 eV, which compares rea- 
sonably well to the experimental value of 1 .52 e V 2 ^. Because 
the HSE06 density functional produces accurate band gaps, 
the accompanying DOS can serve as a plausible representa- 
tion of a photoemission spectrum. Here, the Brillouin zone is 
sampled on a 2 x 2 x 2 T-centered grid and a Gaussian smear- 
ing of 0. 1 eV is applied to orbital occupations to reduce sam- 
pling artifacts in the DOS. For comparison, the HSE06-based 
DOS of crystalline GaAs is computed in VASP with high en- 
ergy resolution using a single unit cell with a 32 x 32 x 32 
sampling of the Brillouin zone and tetrahedral Brillouin zone 
integration. 

III. ELECTRONIC MODELS OF DAMAGE CLUSTERS 

We develop a damage cluster model in three stages of in- 
creasing physical complexity and number of free parameters. 
The goal is to produce a total energy model of a damage clus- 
ter, E(Q), as a function of its total charge, Q, relative to the 
uncharged state, E(Q = 0) = 0. Q is an integer in units of |e|, 
the magnitude of electron charge. Given an electron chemical 
potential, /i, the stable charge state of the cluster, Qs(/x), is 
defined by minimizing the energy over the transfer of charge 
between the crystalline bulk and the cluster, 

E(Q 8 (ji)) + MQs(M) = mm [E{Q) + M Q] . (1) 

Qs(n) is usually positive when fj, is near the valence band 
edge, Ey, and negative when /i is near the conduction band 



edge, Sc- 
Charge transfer energy levels of the damage cluster, Aq,, 
are the values of fj, at which two different charge states, Q and 
Q', have equal energy (Q' > Q by convention). Aq, has a 
simple analytic formula, 

Q _ E(Q)-E(Q') 

Q'-Q ■ () 
The most physically important energy levels are the transi- 
tions between stable charge states, where Q — Qs(Aq, + <5) 

and Q' — Qs(Aq, — 8) for an infinitesimal energy S. Typi- 
cally, the transition between stable charge states is a single- 
electron process with Q' = Q + 1. If the transition is 
a multi-electron process, then all the single-electron Q" to 
Q" + 1 transitions involving metastable charge states between 
Q and Q' are also physically important because they occur 
at a higher rate than the multi-electron process. The simplest 
case is Q' = Q + 2, which is characterized as negative-?/ 
behavior^. Here, the Ag +2 transition between stable states is 

halfway between the Ag +1 and Ag^ energy levels. We only 
consider damage cluster models with single-electron transi- 
tions between stable charge states, which is valid when struc- 
tural relaxation effects are small. 

The first stage of electronic model construction is a basic 
argument for minimal structure. The energy levels of point 
defects can occur anywhere within the band gap between the 
valence band edge and conduction band edge of GaAs, which 
is an energy window of 1.52 eV 2 ^! While computed point 
defect levels in GaAs are observed to have some preference 
for the lower half of the band gap 8 , levels were calculated to 
span the band gap (despite a nominal band gap problem in 
the LDA and PBE functionals), and we assume no bias in the 
larger configuration space of damage clusters. If the distribu- 
tion of levels within the gap is uniform, then it can be charac- 
terized on average by a capacitive total energy model that is a 
quadratic polynomial in charge, 

E{Q) = -noQ + jcQ 2 , (3) 

Aq +1 = mo - ij(Q + |), 

with two free parameters, /io and C, that set an origin and 
energy density for the defect levels. As yet, the model does 
not have any microscopic physical content. 

The second stage of electronic model construction is a 
classical picture that assigns physical meaning to /io and C 
in Eq. We assume that heavily damaged GaAs behaves 
electronically like a featureless metal that is characterized by 
its Fermi level, Ep. This assumption is based on the fact that 
Ga is a metal and As is a metalloid as elemental solids, and 
the forbidden energy gap only emerges from crystalline order. 
The remaining crystalline GaAs is modeled as a featureless di- 
electric medium with a static dielectric constant of e = 12.!^. 
Geometrically, we assume damage clusters to be spheres char- 
acterized by their radius, R. The cluster model is based on the 
energy of a charged metal sphere embedded in a dielectric, 

E{Q) = -E F Q+^ n Q\ (4) 
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in Hartree atomic units. From this model, we hypothesize that 
fio = Ep is a constant for damage clusters in the high-damage 
limit, and C = eR is specific to each cluster and proportional 
to the cube root of its volume. 

The third stage of electronic model construction retains 
the simple geometric picture of the second model and refines 
the electronic behavior of the damage clusters. While suffi- 
ciently high-energy damage to a GaAs crystal should produce 
metallic behavior, low-energy amorphous GaAs is observed to 
retain semiconducting behavior with a reduced band gafPH In 
this case, the common Fermi level for electrons and holes in 
Eq. Q is replaced by valence and conduction band edge en- 
ergies of the semiconducting damage cluster, E a y and E a c, 
with an assumed alignment with respect to the crystalline band 
edges of 



Ey < E a y < Ep < E a c < Ec- 
The updated cluster model is 

E(\Q\) = -E aV \Q\ 



E(-\Q\) = E aC \Q\ 



(5) 
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which continues to assume that all excess charge lies on the 
surface of the sphere as in a metal. While Ep is assumed 
to be universal for a highly damaged region of GaAs, E a y 
and E a c should vary between clusters and approach Ey and 
Eq as damage is annealed and the semiconducting crystalline 
state recovers. 

The assumption of spherical geometry does not reduce the 
generality of the cluster model. We demonstrate this by as- 
suming the damage cluster to be a distribution of point charges 
that have a continuously tunable positive charge, qi. For a set 
of interpoint distances, ry, the minimum Coulomb energy of 
this cluster with a total positive charge of Q is 
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(7) 
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which is a quadratic programming problem. We rescale the 
variables to X{ — q JQ and reduce the problem to 



iQ 2 jgs E^. 



(8) 



1/R 



with a solution that defines an effective spherical radius for an 
arbitrary spatial distribution of point charges. 



IV. RESULTS 

We validate and calibrate the damage cluster model with 
the combined results of MD and DFT simulations. The lim- 
its on time and length scales of direct simulation necessitate 
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FIG. 1. Electronic density of states of crystalline GaAs (gray shaded 
region) compared to amorphous GaAs structures. The structure in 
(a) is formed by annealing at 2500 K. The structure in (b) is the 
first structure with As-As and Ga-Ga bonded minimized by pairwise 
atomic rearrangements and further annealing at 800 K. The Fermi 
level of the amorphous structures is aligned to the center of the crys- 
talline GaAs band gap. 



finite-size scaling assumptions to provide a complete physical 
picture. 

We assume that the two MD-generated bulk amorphous 
GaAs structures are representative of the interior of large dam- 
age clusters that cannot be simulated directly with DFT. Con- 
tingent on this assumption, we verify that heavily damaged 
GaAs behaves electronically like a metal, as observed in pre- 
vious studies 25 . The electronic DOS of the two structures are 
plotted in Fig. [T] against the electronic DOS of semiconduct- 
ing crystalline GaAs. There is little variation between the two 
cases, suggesting an approximately universal electronic DOS 
characteristic of heavily damaged GaAs. Both samples are 
verified to be metallic, but with a strong pseudogap behavior 
that produces a DOS minimum at the Fermi level. 

We estimate the Fermi level of heavily damaged GaAs by 
fitting the models developed in Section [ill] to the vertical de- 
fect levels of the artificial damage clusters described in Sec- 
tion [n] In Fig. [2^, each cluster is independently fit to the 
capacitive model in Eq. (|3j. Each calculated defect level 
A|7q5 is labeled by x, and (1q and C are determined by a 
least-squares minimization of ^{x — x) with 
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(9) 



The optimized values of /xo are approximately normally dis- 
tributed about Ey + 0.46 eV with a standard deviation of 0.07 
eV, as shown in Fig. |2j3. This supports the assumption that the 
electronic properties of these clusters derives from a common 
metallic state with a well-defined Fermi level. 

In Fig. [3^, the clusters are fit to the metal sphere model 
in Eq. (Wl), by least-squares minimization of ~ (x — x) and a 
common value of Ep for all clusters. The optimized value of 
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FIG. 2. Artificial damage clusters fit to a capacitive model of their 
defect levels, AjjjTg'f. (a) The correspondence between x and the fit 
value, x, defined in Eq. |5). (b) The distribution of [1q values from 
the optimized models and a Gaussian fit. 



Ep is 0.46 eV, in agreement with the mean of the /j,q distri- 
bution of the previous fit. With fewer overall free parameters, 
the root-mean-square (RMS) error in x — x grows from 0.31 
to 0.43. The clusters are constructed with a well-defined ge- 
ometric radius which should be comparable to the radius, R, 
contained in the metal sphere model. The geometric radius is 
consistently larger than R, which we attribute to a finite skin 
depth of the metal, i? s kin- A least-squares fit of the electro- 
statically derived radius, R + i? s km, to the geometric radius 
produces a skin depth of 2.1 bohr and an RMS error of 1.3 
bohr, shown in Fig. [3j>. i? s ki n is in reasonable agreement with 
the GaAs skin depth of 1 .6 bohr used to extrapolate the total 
energy of charge defects. To assess finite-size errors of the Ep 
estimate, the metal sphere model is refit to subsets containing 
clusters of fixed geometric radius. No significant dependence 
of Ep on radius is detected in Fig. [3J;. 

We further validate the metal sphere model by visualiz- 
ing the distribution of excess charge within damage clusters 
in Fig. [4] One each of the smallest (8 bohr) and largest (16 
bohr) radius artificial clusters is selected randomly and one 
electron is added to the neutral state. In each case, the defect 
state that the excess electron occupies is delocalized over the 
cluster and has a significant tail extending outside the cluster. 
Because of the large polarizability within the metallic cluster, 
charge fluctuations are heavily screened from the interior and 
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FIG. 3. Artificial damage clusters fit to a metallic sphere model of 
their defect levels, A^q'j. The correspondence (a) between x and 
the fit value, x, defined in Eq. l|9j, and (b) between geometric and 
electrostatic radii, (c) The variation of Ef with cluster size plotted 
against the globally optimized value (dashed line). 
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FIG. 4. Radial distribution of excess charge in two damage clusters 
with one electron added to each. We compare the charge density of 
the defect state (gray shaded region) to the total electronic charge 
rearrangement between the charged and neutral clusters (line) and 
the geometric radius (vertical solid line) to the electrostatic radius 
(vertical dotted line). The electrostatic radius is defined here as the 
minimum radius containing the total excess charge of magnitude e _1 
that remains after screening from bulk GaAs. 



a net screened charge of magnitude eT 1 is confined to within a 
few bohr of the geometric radius of the clusters, in agreement 
with the skin depth fit to the metal sphere model. Charge re- 
arrangement beyond 20 bohr from the center of the clusters 
is omitted from Fig. |4]because it contains finite-size artifacts. 
The LMCC method removes the effect of these artifacts from 
the total energy, but they remain in the electronic charge den- 
sity at the edge of the defect-centered Wigner-Seitz cell. 

The artificial damage clusters and bulk amorphous struc- 
tures have related metallic electronic properties, which is as- 
sumed to derive from sufficient damage to the crystalline 
structure. We assess the degree of damage in these structures 
by the excess cohesive energy per atom of the damaged re- 
gions relative to crystalline GaAs. For amorphous GaAs, typ- 
ical numbers are 0.59 eV/atom for a highly disordered metal- 
lic structure and 0.36 eV/atom for an annealed semiconduct- 
ing structure 25 . The bulk amorphous GaAs structures that 
we have generated have excess energies of 0.61 eV/atom and 
0.58 eV/atom. With only atoms inside the bounding spheres 
counted as damaged, the artificial clusters have an average ex- 
cess energy of 0.74 eV/atom with an RMS deviation of 0.12 
eV/atom between clusters. By this criterion, the clusters are 
more damaged than the bulk amorphous structures, which is 
consistent with the absence of an MD-based annealing step in 
their creation. Physical damage clusters should also have a 
high excess cohesive energy because of kinetic limitations to 
damage annealing. 

There are not enough MD-generated damage clusters to 
produce useful statistics, but we examine them individually to 
identify common behaviors. The details of these clusters are 
summarized in Fig. [5] The vertical defect levels of the clusters 



are fit to the semiconducting sphere model in Eq. ([SJ. There 
is only one significant outlier in this fitting process, Fig. |5}. 
This structure has a large energy window of stability for the 
Q = — 1 charge state instead of the neutral state as predicted 
by the semiconducting sphere model. We interpret this as a 
hole doped into the cluster, resulting from the overall excess of 
5 As atoms. By adjusting for this in the model, the fit improves 
to an RMS error in the energy levels of 0.025 eV. 5 of the 10 
clusters exhibits a band gap, E a c ~ E a y, greater than 0.2 eV, 
which validates the need for a generalization of the metallic 
sphere model. The assumed energy alignment in Eq. <|3j is 
satisfied within the uncertainty of Ep . 

The lifetime of a charged state is likely to be longer than 
the timescale of structural relaxation. Therefore, the adiabatic 
defect levels are more physically relevant than the vertical lev- 
els. However, the relaxation of charge states exacerbates defi- 
ciencies in MD-based annealing of the damage clusters. MD 
simulation does not account for electronic charge fluctuations, 
and only the energy surface of the neutral structure is sam- 
pled. There can exist structural relaxations that significantly 
and irreversibly lower the energy with a large energy barrier 
in the neutral state, but with a small barrier or none at all for 
a charged state. A more physical annealing process would 
enable charge fluctuations that significantly increase the tran- 
sition rate into such lower energy structures. 4 of the 10 MD- 
generated clusters have unstable charge states that cause irre- 
versible structural relaxations. In two cases, Figs. [SJ 7 and[5ji, 
relaxation effects produce a small perturbation in the defect 
levels. In the remaining cases, relaxation effects are large and 
cause many charge states to be unstable at all values of the 
chemical potential. While MD-based annealing does not pre- 
vent the charging instability of the cluster in Fig. |5j or the 
large relaxation effects of the cluster in Fig. [5], it does elim- 
inate the large relaxation effects of the cluster in Fig. |5^ by 
transforming it into the cluster in Fig. [5}i. This evidence is 
not conclusive, but it demonstrates that MD-based annealing 
is capable of reducing relaxation effects but incapable of reli- 
ably relaxing unstable charge states. Relaxation effects were 
not studied for the artificial damage clusters. We conjecture 
that more thorough annealing would significantly reduce re- 
laxation effects on average and that our study of vertical de- 
fect levels is representative of the adiabatic defect levels of 
clusters that have undergone such an annealing process. 



V. DISCUSSION 

The results of the previous section clarify that direct MD- 
plus-DFT simulation of damage clusters in GaAs are limited 
by both the size of clusters studied in DFT and the timescale of 
clusters annealed in MD combined with a lack of charge fluc- 
tuations in MD. However, the electrical charging of even very 
small clusters are accurately fit to a classical model, and we 
expect the model to continue to work for larger, "more classi- 
cal" damage clusters. The model contains a geometric factor 
that can be estimated without quantum mechanical simulation 
and amorphous semiconductor band edges that asymptote to 
a universal Fermi level in the high-damage limit. What re- 
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FIG. 5. 10 GaAs damage clusters extracted from MD simulations, with (h-j) further annealed after extraction ((h) is the annealed version of 
(g)). Structures are visualized in QuteMol 29 , with Ga and As as black and white spheres. Only atoms that deviate from their ideal crystalline 
positions by more than 5% of the lattice constant are shown. Total number of defective atoms in the cluster are given, with +/— denoting the 
deviation from the crystal stoichiometry. The fit to the semiconducting sphere model in Eq. j6j as well as the RMS error in Aq +1 is given 
for each cluster. Transitions between stable charge states, Ag, , are shown along with the Q value of the energy intervals for the model (M), 
vertical (V), and adiabatic (A) calculations of cluster energies. Q' = Q + l for all model and vertical charge states, and only the Q — charge 
state is labeled. 



mains unknown are the initial population and size-distribution 
of damage clusters as a function of incident ion energy, the 
long-time-dependence of cluster annealing, and the change in 
electronic properties as damage within the cluster is repaired 
by annealing. Characterization of the initial distribution of 
damage clusters will be reported in a separate papeJ^l The 
remaining unknowns require further work to completely clar- 
ify. 

It is worthwhile to compare our partially constructed 
model of damage clusters with the existing radiation damage 
model of point defects. A key difference is that the number 
of defect levels in a dilute distribution of point defects is pro- 
portional to the number of damaged atoms, while the number 
of levels in a damage cluster is dependent on the volume of 
the cluster and insensitive to the number of damaged atoms, 
as apparent in Fig. [3] Superficially, this is a significant dif- 
ference between models. However, the point defect model 
adjusts the energy levels of the point defects with the local 
electrostatic potential, which captures the same basic effect 
as the electrostatic charging energy of the models developed 



in Section III for a dense distribution of charged points de- 
fects. While damage clusters cannot be structurally decom- 
posed into well-defined point defects in general, it may be 



possible to produce the same set of defect levels with a care- 
fully chosen high-density distribution of point defects. With 
such an approximation scheme, the electronic behavior of a 
multiple-point-defect model is able to mimic cluster behavior. 

Because a radiation damage model based on only point 
defects is capable of approximating the electronic properties 
of damage clusters, fitting such a model to experiments may 
prescribe an excess population of point defects to produce the 
same electronic behavior as the missing clusters. Whether or 
not such an effective model is accurate depends on the dif- 
ference in time-dependence between a dense distribution of 
point defects and a damage cluster. The rate of binary chem- 
ical reactions between point defects will increase with den- 
sity, which might reduce the lifetime of a dense distribution 
of point defects relative to a damage cluster. The long-time 
evolution of a damage cluster is not yet available for compar- 
ison. Also, the density of a distribution of point defects can 
be dynamically reduced either by diffusion or by drift result- 
ing from Coulomb repulsion between charged point defects. 
These are dynamical degrees of freedom that a damage clus- 
ter is not likely to have, assuming some degree of cohesion 
within the cluster and negligible mobility of the cluster as a 
whole. 



8 



We propose a simple rationalization of the time evolution 
of a damage cluster, which directly relates to the parameters of 
the semiconducting sphere model in Eq. d3]l. First, we assume 
that the amorphous GaAs band edges are directly correlated 
to the average excess cohesive energy per atom over the crys- 
talline state. Next, we assume that the effect of annealing can 
be decomposed into a reduction of the excess cohesive en- 
ergy plus a surface recrystallization that is nucleated by the 
crystalline region outside the cluster. These effects are a ra- 
tionale for the distinct time-dependence of E a v and E a c and 
of R. The rate of excess cohesive energy reduction and the 
speed of the crystallization front can be studied separately in 
idealized simulations of amorphous GaAs or of an interface 
between crystalline and amorphous GaAs. In addition to the 
evolution of the damage cluster, there also may be dynamic 
interactions between clusters and point defects. For example, 
damage clusters might act as a source or sink for point defects. 
The exploration of these ideas is left to future work. 



VI. CONCLUSIONS 

A complete radiation damage model in GaAs must ac- 
count for the initial conditions, electronic properties, and 
time-dependence of all relevant material defects following a 
radiation event. It has been assumed previously that all rel- 
evant defects are simple, enumerable point defects, but we 
have established the commensurate relevance of larger dam- 
age clusters. The defect energy levels of all damage clusters 
cannot be exhaustively computed, but they are well described 
by a classical total energy model fit to a statistical sampling 
of damage clusters. This result highlights the importance of 
sampling a statistical distribution when studying physical phe- 
nomena involving disorder at the atomistic scale. A com- 
plementary report on the initial distribution of damage clus- 
ters observed in MD simulations of collision cascades is in 
preparation 17 . We have outlined a possible model of the long- 
time evolution of damage clusters, but validation and parame- 
terization of the model remain as future work. This remaining 
problem is difficult because of the long timescales involved 
and the possibility that electronic charge fluctuations are an 
important relaxation mechanism not present in tools such as 
MD that are capable of directly accessing the necessary time 
scales. 
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